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ABSTRACT 

In a previous paper, we modeled the oscillations of a thermally-supported (Bonnor-Ebert) sphere 
as non-radial, linear perturbations following a standard analysis developed for stellar pulsations. The 
predicted column density variations and molecular spectral line profiles are similar to those observed in 
the Bok globule B68 suggesting that the motions in some starless cores may be oscillating perturbations 
on a thermally supported equilibrium structure. However, the linear analysis is unable to address 
several questions, among them the stability, and lifetime of the perturbations. In this paper we 
simulate the oscillations using a three-dimensional numerical hydrodynamic code. We find that the 
oscillations are damped predominantly by non-linear mode-coupling, and the damping time scale 
is typically many oscillation periods, corresponding to a few million years, and persisting over the 
inferred lifetime of gobules. 

Subject headings: stars: formation, hydrodynamics, radiation mechanisms: non-thermal, line: profiles, 
ISM: globules, ISM: clouds 



1. INTRODUCTION 

Recent observations suggest that the internal structure 
of some of the small isolated molecular c louds known as 
starless cores or Bok globules (|Bok 1119481 ) is well approx- 
imated by a balance of thermal and gravitational forces 
on which is superposed a pattern of non-radial oscilla- 
tions (Lada et al. 2003). Previously, the internal oscilla- 
tions of these sta rless cores have b een noted in simula- 
tions of coll apse (|Hennebelle II2003I ). magnetized clouds 
! (iGalh 1120051) . modeled as one-dimensional radial oscilla- 
tions by means of a numerical simulation that includes 
radiative energy losses (Keto & Field 2005) and as three- 
dimensional no n-radial perturba tions assuming isother- 
mal conditions (jKetolitl^lllQQ^ . Analysis of these mod- 
\ els and comparison of simulated molecular spectral line 
profiles with those observed in the dark clouds L1544 
■ (jCaselh et al| [2002) and B68 ()Lada e t al. 2003) indicate 
' that the motions in some dark clouds are consistent with 
large-scale hydrodynamic oscillations, i.e., sound waves. 
Nevertheless, a number of questions remain as yet unan- 
swered. 

Is the underlying assumption of the oscillations as per- 
turbations valid for the amplitudes required to produce 
the observed velocities? Based on the observations of 
B68, the perturbation amplitudes in our previous three- 
dimensional model were ~ 25%. Does the linear model 
correctly predict the velocity field for such large ampli- 
tudes or would significant differences appear in the non- 
linear regime? 

What is the lifetime of the perturbations and thus the 
velocities within the clouds? For hydrodynamic oscil- 
lations to remain a viable explanation of the velocity 
field observed in B68 they must have lifetimes in excess 
of a crossing time. Supersonic oscillations are strongly 
damped via shocks. However, in the subsonic regime 
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the dominant damping mechanism is less clear. In the 
one-dimensional model (|Keto fc Fieldl[2005[ ) . the dissipa- 
tion of the oscillations occurs via radiative losses with an 
estimated timesc ale of 10 Myr. I n the non-radial per- 
turbative model (|Keto et al.ll2006l ) damping is explicitly 
ignored. For mode amplitudes larger than linear per- 
turbations, non- linear mode-mode coupling provides an 
additional potentially significant damping mechanism. 

Are the oscillations stable or growing? Are there pre- 
ferred modes? In general we expect the oscillations to 
damp over time, and we expect that the longest-lived 
modes of oscillation would be the long wavelength modes. 
Neither the one-dimensional non-linear analysis nor the 
three-dimensional linear analysis is able to adequately 
address this question. 

In this paper we address these questions by modeling 
the oscillations of dark clouds with a three-dimensional, 
non-linear, numerical hydrodynamic simulation. We find 
that the dominant damping mechanism is indeed non- 
linear mode coupling, the longest wavelength mode life- 
times are on the order of a few Myr and produce velocity 
and density fields qualitatively similar to those observed. 
Section [2] summarizes the numerical methods employed, 
sections [3] and HI discuss the viability of quadrupole and 
dipole oscillation models, and the conclusions are con- 
tained in section [5l 

2. COMPUTATIONAL MODELING 

2.1. Hydrodynamic Evolution 

Because the thermal gas heating and cooling time, 
via co Uisional coupling to dust at high densities 
(jBurke & HoUcnbach 1983) and molecular line radiation 
at low densities (|Goldsmith Il200lh . is short in compar- 
ison to the typical oscillation period (on the order of 
10^ yr), the oscillations of dark-starless cores are well 
modeled with an isothermal equation of state. In par- 
ticular, we use a barotropic equation of state with an 
adiabatic index of unity, in which case the isothermal 
evolution is also adiabatic. Since the unbounded isother- 
mal gas sphere is unstable, we truncate the solution at a 
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Fig. 1. — The density (solid) and pressure (dashed) profiles of 
the unperturbed isothermal cloud (Bonnor-Ebert sphere). Near 
the surface of the cloud {r/R = 1) the density rapidly decreases. 
However, the transition region in the equation of state softens this 
drop over many grid zones to avoid difficulties in resolving the 
pressure gradients. 
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Fig. 2. — The dimensionless amplitude of the first quadrupole 
p-mode, the (1,2,2) mode. The dashed line is a quadratic fit to the 
logarithm of the amplitude, with an e-folding time of roughly 2 Myr 
at t = 0. The squares indicate the times for which the frames in 
Figures |3] were computed. The dotted and long-dashed lines show 
the only other two significant long-wavelength oscillations, listed 
by their nlm. 



given density contrast {p{R) / pc), producing the standard 
Bonnor-Ebert sphere solution. This is done by inducing 



a phase change in the equation of state, i.e. 
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where the intermediate state, with 6 — 0.8 and e = 
log((5)/ log((5T/Text) ^ 1 chosen to make P{p) contin- 
uous, is introduced to avoid numerical artifacts at the 
surface, and Text 3> T. The resulting radial density and 
pressure profiles are shown in Figure [T] Note that there 
is a transition region, associated with the intermediate 
state in the equation of state, which distributes the sud- 
den drop in density over many grid zones. Despite this, 
the pressure is continuous (necessarily) and nearly con- 
stant outside the surface of the Bonnor-Ebert sphere. We 
find that generally our results are insensitive to the par- 
ticular equation of state of the exterior hot gas, given 
that it is sufficiently dilute. This is appropriate for glob- 
ules like B68 which are known to be surrounded by hot, 
low density gas. 

We employ a three-dimensional, self-gravitating hydro- 
dynamics code to follow the long term mode evolution. 
Details regarding the numerical algorithm and specific 
validation for oscillating g as spheres (in that case a whit e 
dwarf) can be found in iBroderick fc Rathore I (|2006f l. 
and thus will only be briefly summarized here. 

The code is a second order accurate (in space and time) 
Eulerian finite-difference code, and has been demon- 
strated to have a low diffusivity. Because the equation of 
state is barotropic, we use the gradient of the enthalpy 
instead of the pressure in the Euler equation since this 
provi des better stability in the unperturbed configura- 
tion fjBroderick fc Rathore l[2006f l. The Poisson equation 
is solved via spectral methods, with boundary conditions 
set by a multipole expansion of the matter on the com- 
putational domain. 

In order to assess the suitability of the code for ap- 
plication to isothermal gas spheres we tested the stabil- 
ity to monopolar perturbations as a function of central- 
to-surface density contrast, finding instability for values 
above 13.5. This is in approximate agreement with the 
rediction of instability fo r density constrasts above 14.3 
Foster fc Chevalier I [l993l) . The deviation from 14.3 is 
likely due to the effective atmosphere produced by the 
exterior hot matter, and in particular the material in the 
transition region of the equation of state. 

Henceforth, the specific physical parameters of the un- 
perturbed model that we use, chosen to be roughly appli- 
cable to B68, are pc = 1.185 x 10"i**gcm-^ p{R)/pc = 
12.5, and T = lOK. In all cases the simulations were 
performed at resolutions of 64'^, 128'^ and 256'^. Com- 
parison of the latter two indicated that the numerical 
solution had converged by these resolutions. The results 
presented in the following sections are from the simula- 
tions at 256^. 

The initial perturbed states are computed in the linear 
approximation using the standard form alism of adiabatic 
(which in this case is also iso thermal, iKeto et aTl 120061 ) 
stellar oscillations (|Coxlll980f) . The initial conditions of 
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Fig. 3. — The equatorial gas density, velocity, and CS(2-1) line profiles of the quadrupolar mode after the first half-period, (left), and 
the second full period, (right). These correspond to the times marked by black squares in Figure [2] Despite having oscillated for nearly 
5 X 10^ yr, the quadrupolar structure is still clearly evident in the density & velocity fields as well as in the spectral line profiles. Spectral 
line profiles colored blue (red) indicate approaching (retreating) gas on the near side relative to the observer. 



the perturbed gas sphere are then given by the sum of 
the equilibrium state and the hnear perturbation. We 
employ the same definition for the d imensionless mod e 
amplitude as given in the Appendix of lKeto et an(|2006f ). 

Throughout the evolution, the amplitude of a given 
oscillation mode, denoted by its radial (n) and angular 
(/ & m) quantum numbers, may be estimated by the 
explicit integral: 

Anim = J d^a; p V • ^li^ , (2) 

where v is the gas velocity, £„;„j is the mode displace- 
ment eigenfunction, ujnim is the mode frequency and the 
dagger denotes hermitian conjugation.^ 

2.2. Molecular Line Simulations 

Assuming that the gas contains the usual interstel- 
lar molecules, we simulate molecular line observations 
of the hydrodynamic model by calculating the non-LTE 
(non-thermodynamic equilibrium) emission of the CS(2- 
1) transition using our three-di mensional ALI (acceler- 
ated A iteration) code (jKeto et al.l I2004J ). The densi- 
ties, velocities, and constant temperature are taken from 
the hydrodynamic simulation. The line broadening is 
assumed purely thermal. 

The abundance of CS is based on a, simp lified model 
of CO chemistry from lKeto fc Casellil (|2008l ). We briefly 
describe the model here. The abundance of CO in the gas 
is decreased by two effects, by photodissociation at the 

^ For stellar oscillations the mode amplitudes may be determined 
from the density perturbation alone using the velocity potential. 
However, when there is a non-vanishing surface density, as is the 
case for the Bonnor-Ebert sphere, this is not possible and the ve- 
locity field must be used. 



edge of the cloud and by freezing onto dust grains in the 
high density center. We calculate the reduction in the 
CO abundance due to these two affects and assume that 
the s ame reduction affe cts the CS molecule. We assume, 
as in lKeto et al.l (|2006f ). an undepl eted abundance of C S 
of 6 X 10"^. This is consistent with lTafalla et al. I (|2002D . 
who cite a range in values derived from their observations 

o f 3-9 X IQ-^ 

iTielens fc HollenbachI (|1985f ) suggest that at the edge 
of the cloud where photodissociation is important, the 
dominant cycle for the formation and destruction of CO 
is 

C+ ^ CH2 ^ CO ^ C ^ C+ . (3) 

The time scale for the formation of CH2 by radiative as- 
sociation is longer than for the formation of CO from CH2 
and O. Therefore we assume that the rate of formation of 
CO is given by that of CII2. The CO cycle may then be 
described by three rate equations for the creation of 0+, 
C, and CO plus the conservation condition that the to- 
tal amount of carbon is contained in these three species. 
The rate equations for photodissociation depend on the 
mean visual extinction at each point, which we calculate 
using the machinery of the three-dimensional radiative 
transfer code. 

The distribution of CO between the gas and grain 
phases is controlled by the rate of depletion from the 
gas phase by freezing onto dust grains and the rate of 
desorption off dust grains by cosmic rays. The steady 
state abundance of CO in the gas phase is then given by 
the ratio of the depletion time to the sum of the deple- 
tion and desorption times. Because the depletion time 
scales inversely with the collision rate while the desorp- 
tion time scales inversely with the cosmic ray ionization 
rate, the depletion of CO is dependent only on the gas 
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Fig. 4. — A map of the CS(2-1) line profile overlaid upon the 
gas column density for the quadrupolar oscillation after nearly 
5 X 10^ yr. This is viewed from A = 30°, 6» = 30°, for which 
the pattern of spec tral line profiles is similar to the Bok globule 
B68 (cf. Figure 7 of lKeto et al|[2006h . The gas density is shown as 
a filled contour map, with contours equally spaced in 10% incre- 
ments of the maximum value and the 50% contour is marked by the 
thick black line. The CS(2-1) line profiles are associated with lines 
of sight centered in each subplot. Blue (red) lines correspond to 
approaching (retreating) gas on the nearside of the cloud relative 
to the observer. Note the prominent double reversal of the line-of- 
sight motions and the substantial distortion of the cloud shape, as 
seen in B68 (and in the same relative orientation). 

density because the cosmic ray ionization rate is the same 
throughout the cloud. 

3. QUADRUPOLE PERTURBATIONS 

The primary mode discussed in lKeto et al.1 (|2006f ) was 
the first quadrupole p-mode, i.e., the n=l, l = m = 2 
oscillation mode. In order to reproduce the observed ve- 
locities it was necessary to consider a dimensionless am- 
plitude of 0.25, which produced near-vanishing densities 
in some portions of the cloud, and is thus clearly well 
into the non-linear regime. 

The evolution of the dimensionless amplitude of the 
(1,2,2) mode, as defined by equation ([2|), is shown as 
a function of time in Figure [2] The decay is well ap- 
proximated by an exponential. However, the strongly 
non-linear nature of the decay is evident in the changing 
decay constant (resulting in a noticeable bowing of the 
fitted line in Figure [2|) , which yields an initial e- folding 
time of approximately 2 Myr, corresponding to roughly 
8 oscillation periods (16 pattern periods). Since the in- 
trinsic dissip ation in the code has previou sly been shown 
to be small (jBroderick fc Rathore Il20()6l ). the decay in 
the amplitude is dominated by non-linear mode interac- 
tions, resulting in the transfer of energy from large-scale, 
low-order modes to smaller-scale, higher-order modes ^. 
However, mass loss driven by large velocity perturba- 

^ Generally the amplitude evolution may be expected to depend 
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Fig. 5. — The dimensionless amplitude of the first dipole p-mode, 
the (1,1,1) mode. The dashed line is a quadratic fit to the logarithm 
of the amplitude, with an e-folding time of roughly 2 Myr at t = 0. 
The squares indicate the times for which the frames in Figure \E\ 
were computed. The long-dashed line shows the (0, 2, 2) mode, the 
only other significant oscillation. 



tions and mode coupling in the artificial atmosphere (see 
Figure [1]) also contribute to the decay. It is not clear 
if these processes are artifacts of our numerical simula- 
tion and the manner in which we initially prepare the 
cloud or if they may be found in nature as well. Thus 
our estimate of the lifetime should be considered a lower 
limit. 

The only two oscillations of comparable wavelength 
(less than approximately R/5) that reach significant am- 
plitudes are the (0,2,2) and (0,4,4) modes, whose cou- 
plings are likely enhanced by the Cartesian geometry of 
the underlying computational grid. In particular, there 
is no evidence for an inverse-cascade, implying that the 
gas cloud will not fragment. Indeed, this can be explicitly 
seen in the density and velocity structures in the equa- 
torial plane (see Figure [3]) in which the mode energy is 
primarily dissipated into small scale motions. Thus even 
vigorous oscillations (at least up to amplitudes of 25%) 
do not appear to be sufficient to split the cloud into a 
binary. 

Despite the decay in amplitude, the quadrupole mode 
structure is still apparent in the density distribution of 
the cloud after a number of oscillation periods. This is 
also apparent in the substantially distorted H2 column 
densities, shown in Figured] for a viewing angle that re- 
produces the generic features of B68 {d and A correspond- 
ing to the lattitude and longitude, respectively, of the line 
of sight and oriented such that A = 0°/90° correspond- 
ing to a view from the bottom/left of the page in Figure 

upon the initial mode spectrum. We investigated the significance 
of this dependence by introducing higher order modes with ran- 
domly chosen amplitudes and phases. We found that as long as 
the quadrupole mode remained dominant (in terms of A„;„i) the 
results were unchanged. 
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Fig. 6. — The equatorial gas density, velocity, and CS(2-1) line profiles of the dipole mode after the first half-period, (a), and the second 
full period, (b). These correspond to the times marked by black squares in Figure [5] Despite having oscillated for more than 5 X 10^ yr, 
the dipolar structure is still clearly evident in the density & velocity fields as well as in the spectral line profiles. Spectral line profiles 
colored blue (red) indicate approaching (retreating) gas on the near side relative to the observer. Note that the spectra along the bottom 
would seem to imply expansion while those along the top would imply collapse. Furthermore, spectra from the side suggests rotation. Thus 
inferring the global state of the cloud from line profiles may be complicated by the presence of large-scale oscillations. 



[31). Even after many oscillations the column density maps 
are generally qualitatively similar to observations of dark 
cores which often s how clouds with an aspect ratio 
of approximately 2:1 (I Ward-Thompson. Motte. &: Andrei 
ll999HBacmann et al. HIoM k 

Similarly, the complex velocity patterns, as inferred 
from the shape of the CS(2-1) line profiles, are present 
after multiple oscillations. To the left and below the 
plots of the density and velocity structure in Figure [3l 
the CS(2-1) line profiles are shown for a number of lines 
of sight. Comparison of the early time and late time line 
profiles suggests that the small scale motions, the result 
of damping of the initial oscillation, do not significantly 
alter the line shapes or, as a direct consequence, the in- 
ferred line-of-sight velocity distribution throughout the 
cloud. Due to the intrinsically complex nature of the 
oscillation velocity field, the inferred velocity map can 
mimic expansion, contraction, rotation, or a combination 
of these. For example, the velocity patterns inferred by 
the line shapes shown in Figure [3] appear to correspond 
to rotation, while if these were rotated by 45° (diago- 
nal lines of sight) they would appear to correspond to 
expansion or contraction. 

The structure of the inferred velocity map can become 
substantially more complex as we approach the pole (the 
polar angle 9 approaches 0°). As seen in Figured! for 



9 ~ 30° and azimuthal angle A ~ 30° a quadrupolar 
oscillation can reproduce the general features of B68. 
This includes a double reversal of the line-of-sight ve- 
locity across the cloud (represented by the change from 
red to blue to red in the line profiles) and the signifi- 
cantly distorted gas column density. In fact, for the view 
shown, the relative orientation of the column density dis- 
tortion and the velocity double reversal is similar to that 
observed in B68! This continues to be the case even after 
5 X 10^ yr. 

As anticipated in iKeto et all ()2006l ). the qualitative 
features of the linear column density and velocity maps 
are still present in the non-linear computation. In partic- 
ular, both the anisotropic nature of the column densities 
and the complex velocity structure persist for a number 
of oscillation periods. Nevertheless, the non-linear decay 
of the amplitude, and subsequent small scale motions, 
does become significant on timescales of a Myr as evi- 
denced by the decay of the mode amplitude in Figure |2l 
As a consequence, non-linear mode damping appears to 
be the dominant constraint upon whether such an oscil- 
lation can be clearly observed in velocity maps. 

4. DIPOLE PERTURBATIONS 

Alth ough not explicitly considered in iKeto et al] 
(|2006l ). dipole oscillations are a possible alternative in- 
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Fig. 7. — A map of the CS(2-1) line profile overlaid upon the gas 
column density for the dipolar oscillation after more than 5 X 10''' yr. 
This is viewed from the same angle as shown in Figure]?] (A = 30° , 
d = 30°). The gas density is shown as a filled contour map, with 
contours equally spaced in 10% increments of the maximum value 
and the 50% contour is marked by the thick black line. The CS(2- 
1) line profiles are associated with lines of sight centered in each 
subplot. Blue (red) lines correspond to approaching (retreating) 
gas on the nearside of the cloud relative to the observer. Note 
the prominent single reversal, similar to that associated with cloud 
rotation. For the dipole mode the distortion in the column density 
map is substantially less than that observed for the quadrupolar 
mode. 

terpretation to rotation, and perhaps also capable of ex- 
plaining the velocity structure of B68. For this reason 
we considered the first dipole p-mode, i.e., the n ^ I = 
m = 1 oscillation mode, as well. 

As a consequence of its longer wavelength, the ampli- 
tude of the dipole mode must be larger than that of the 
quadrupole oscillation mode to produce similar veloci- 
ties. Hence, we chose a dimensionless amplitude of 0.5. 
Similar to the quadrupole case, this is also likely to be 
strongly non-linear. 

As in the quadrupole case, the dipole oscillation de- 
cays primarily into small scale motions. In both the 
quadrupole and dipole case the cloud sheds a few per- 
cent of its mass from the outer layers, contributing to the 
decay of the oscillation. In the dipole case this also pro- 
duces a notable shift in the center of mass of the cloud as 
a result of the asymmetric mass loss through the bound- 
ary of the computational domain. 

As with the quadrupolar oscillations the gas column 
densities are initially strongly distorted. However, in 
contrast to the quadrupole mode, after the first oscil- 
lation period (approximately 1.5 x 10^ yr) they have be- 
come roughly circular (see, e.g.. Figure [7]). This suggests 
that large amplitude dipole modes alone are insufhcient 
to generate the observed anisotropy in the shapes of dark 
cores. 

Figure [5] shows the evolution of the dimensionless am- 



plitude for the first dipole p-mode. Again the evolving 
decay constant implies that this is a strongly non-linear 
process. However, in this case no oscillations with similar 
amplitudes develop within the duration of the simulation 
(the largest being the (0, 2, 2) mode, which again is likely 
due to the Cartesian nature of the grid) . The mode life- 
time is on the order of 2 Myr, considerably longer than 
the mode oscillation period. As before, there is no indi- 
cation of cloud fragmentation. 

Nevertheless, the velocity maps inferred from the CS 
(2-1) line shapes, shown in Figures [6l retain their dipo- 
lar structure, even after many oscillations. Note that this 
implies that inferring the global dynamical state of the 
cloud from line profiles alone is significantly complicated 
by the possible existence of large-scale oscillations. As 
shown in Figures [6] & El for a variety of viewing angles, 
the motion appears to mimic rotation, expansion or con- 
traction, all while in the same dynamical state. This is 
still true after the column density maps become nearly 
spherical. 

5. CONCLUSIONS 

Despite their strongly non-linear nature, large- 
amplitude oscillations of pressure-supported molecular 
clouds reproduce many of the qualitative features of ob- 
servations of starless cores. These include the complex 
velocity structures and highly distorted column densi- 
ties found in B68 and L1544, as well as mimicking ro- 
tation, expansion, and contraction depending upon the 
mode structure and viewing angle. This may complicate 
efforts to discern the global dynamical state of these ob- 
jects from observations of molecular line shapes alone. 

However, it may be possible to distinguish oscillations 
from collapse by comparing the spectral line profiles and 
velocities of volatile molecules such as N2H-I- and NHS 
against those of more refractory species such as CO, 
HCO-H, or CS (|Keto fc Fieldl[20?)5l ). In a simple breath- 
ing mode oscillation, the gas velocities are highest to- 
wards the outer regions of the cloud where the abun- 
dance of the carbon species is not depleted by freeze-out. 
In contrast, in the gravitational collapse of unstable BE 
spheres, the gas velocities are highest in the center of the 
cloud where the emission from the non-depleting nitro- 
gen molecules is highest. 

We note that we have restricted our attention to hy- 
drodynamic oscillations of pressure-bounded, thermally- 
supported spheres. This description seems appropriate 
for globules such as B68 that are surrounded by hotter 
gas and have density profiles that closely match those of 
Bonnor-Ebert spheres. Whether the Bonnor-Ebert de- 
scription is appropriate for embedded cores remains to 
be determined. For example, if magnetic stresses con- 
tribute significantly to cloud support, then additional 
classes of osci llations exist, associated with t he magnetic 
field (see.e.g.. iHennebelle 1120031 : [C^aliTI 120051 ). 

The lifetimes of hydrodynamic oscillations of starless 
cores like B68 are sufficiently long (few 10^ yr) to be a 
significant fraction of the expected lifetime of a molecular 
cloud (10^ yr). The oscillation lifetimes are dominantly 
limited by non-linear mode-mode coupling. Nevertheless, 
large amplitude oscillations persist for many periods (and 
thus many sound crossing times), and are therefore not 
expected to be particularly rare. There are no signs of 
cloud fragmentation, and thus amplitudes comparable to 
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those discussed here are insufficient to produce binaries 
without the inclusion of additional physical processes. 
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